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Abstract 

We evaluate the Brownian motion of a massive particle ("black hole") at the center of a galaxy using N- 
body simulations. Our galaxy models have power-law central density cusps like those observed at the centers 
of elliptical galaxies. The simulations show that the black hole achieves a steady-state kinetic energy that is 
substantially different than would be predicted based on the properties of the galaxy model in the absence of 
the black hole. The reason appears to be that the black hole responds to stars whose velocities have themselves 
been raised by the presence of the black hole. Over a wide range of density slopes and black hole masses, the 
black hole's mean kinetic energy is equal to what would be predicted under the assumption that it is in energy 
equipartition with stars lying within a distance ~ r/,/2 from it, where r/, is the black hole's influence radius. 

The dependence of the Brownian velocity on black hole mass is approximately (V^) °^ MgU^^ with y the 
power-law index of the stellar density profile, p r^^. This is less steep than the M^^ dependence predicted 
in a model where the effect of the black hole on the stellar velocities is ignored. The influence of a stellar mass 
spectrum on the black hole's Brownian motion is also evaluated and found to be consistent with predictions 
from Chandrasekhar's theory. We use these results to derive a probability function for the mass of the Milky 
Way black hole based on a measurement of its proper motion velocity. Interesting constraints on Mbh will 
require a velocity resolution exceeding 0.5 km s^'. 
Subject headings: stellar dynamics, galaxies; nuclei, black holes 



1. INTRODUCTION 

A massive black hole at the center of a galaxy undergoes a 
random walk in momentum space as its motion is perturbed 
by gravitational encounters with nearby stars. The expected 
amplitude of this "gravitational Brownian motion" is 

x( ° . y ,1) 

VlOOkms^V 

where M and m are the mass of the black hole and a typical 
star respectively and a is the ID stellar velocity dispersion; 
the brackets denote a time average. Equation Q follows from 
assuming that the average kinetic energy of the black hole is 
equal to that of the stars. Gravitational Brownian motion is 
potentially interesting for a number of reaso ns; as a means t o 
constrain black hole masses (e.g. iBacker & Sramekl J1999I) ; 
iReid et al." (199?, "2003")); as a mechanism for displacing black 
holes from their otherwise central locations in galaxies (e.g. 
iBahcall & W olf ( 1976)); and as a possible source of enhance- 
ment in the rate of supply of stars to the black hole (e.g. lYound 
©80)). 

The approximate validity of equation ([Q 
has been conf irmed i n a number of numer- 
ical studies (iMilled Il992t i Taga & Ive 1998; 
i Milosavlievic & Merrittll20 0lUCh atteriee. Hernauist. & Loeb 
l2002albt Dorba nd. Hemsendorf. & MerrittI 



'2003' [Chatteriee. Hernauist. & Loebl '2003t 

■Makino & Funato 2004. With two excep- 

tions however jMilosaylievic & Merritt 2000 
iDorband. Hemsendorf. & MOTitI |2003Ik these studies 
have been based on galaxy models with large, low-density 
cores, very different from the dense, p ~ r^^ nuclei observed 
at the centers of the galaxies known to harbor supermassive 
black holes. Furthermore the black hole particle is often 
introduced into the simulations in a non-self-consistent way; 
first a black-hole-free model is constructed, then the black 
hole is added, which causes the galaxy model to evolve away 
from its initial state as the stellar motions respond to the sud- 
denly dee pened potential (e.g. Chatteriee. Hernauist. & Loe^ 
J2002al b')'). These practices complicate the interpretation 
of the numerical simulations and make it difficult to derive 
predictions about the expected behavior of black holes in real 
nuclei. For example, the quantity a that appears in equation 
Q is well defined at the center of a low-density galaxy 
containing no black hole, but in a real galaxy, o is a strong 
function of radius, due both to the inhomogeneity of the 
galaxy and to the presence of the black hole itself. 

These considerations motivated us to undertake a new series 
of numerical studies of gravitational Brownian motion. Our 
primary goal was to evaluate the validity of the equipartition 
assumption for black holes at the centers of galaxy models 
with realistically high central densities. We carried out a set 
of A^-b odv integrations similar in character to those described 
bv lborband. Hemsendorf. & MerrittI (i2003l) . but with a wider 
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range of galaxy models and, in some cases, with a spectrum 
of stellar masses. All of the galaxy models have a power-law 
dependence of stellar density on radius near the center, and 
the initial stellar velocities were generated from a distribu- 
tion function that accounts for the gravitational force from the 
"black hole"; hence, the models are in a precisely equilbrium 
state at time zero. 

Our results can be summarized in a surprisingly simple 
way: the black hole reaches a state of energy equipartition 
with the stars lying at distances ^ r/,/2 from it, where r/, is 
the black hole's gravitational influence radius (a precise defi- 
nition of r/i is given below). By definition, stars within r/, are 
moving largely in response to the gravitational force from the 
black hole; hence, the black hole's Brownian motion is deter- 
mined by the velocities of stars which themselves have been 
raised by the presence of the black hole! One consequence is 
that the black hole's rms velocity does not drop as steeply with 
black hole mass as the M^'/^ dependence predicted by equa- 
tion Q. Another is that (V^) can be substantially different - 
higher or lower - than the value predicted by equation Q if o 
is measured outside of the black hole's influence radius. 

The properties of our galaxy models and the details of the 
A^-body integrations are presented in §2 and §3. The predic- 
tions of Chandrasekhar's theory of stellar encounters are re- 
viewed in §4; as we point out there, Chandrasekhar's theory is 
essentially local in character and does not make useful predic- 
tions about the expected value of (V^) in realistic, inhomoge- 
neous galaxies. Nevertheless the theory does say something 
definite about the relation between (V^) in the single- and 
multi-mass cases. Results from the A^-body integrations are 
presented in §5 (single stellar mass) and §6 (mass spectrum). 
In §7, we use our results to predict the expected amplitude of 
the Brownian motion for the Milky Way black hole, and show 
how a measurement of its velocity, or determination of an up- 
per limit, can be converted into a probability function for its 
mass. 

2. GALAXY MODEL 

Supermassive black holes are observed at the centers of 
early-type galaxies and the bulges of spiral galaxies. The 
luminosity profiles of these systems are well represented as 
power laws in the radius at distances from the black hole less 
than ^ r/,, the break radius. (So-called "core" galaxies also 
have luminosity profiles that are power laws in the space den- 
sity, but with indices less than one.) We adopted Dehnen's 
(1993) density profile for our galaxy models: 

(3-Y)Mgai a 
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(2) 



where p (r) is the stellar mass density, Mgai is the total mass in 
stars, a is the scale length (roughly speaking, the break radius) 
and y is the logarithmic slope of the central density cusp. 

To this galaxy model was added a central point of mass M 
representing the black hole. For our purposes, it was crucial 
that the initial model be in a steady state, so that any time de- 
pendence could be attributed to perturbations resulting from 
the finite-A^ realization of the model, and not to ill-defined 
departures from equilibrium. To achieve this, the initial ve- 
locities of the stars were generated from the unique, isotropic 
phase-space distribution function that reproduces the Dehnen 
density law in the smoot h combined pot ential of the stars and 
the black hole particle (ITremaine et al.i ri994). This distribu- 
tion function is non-negative for y > 0.5 in the presence of 
a central point mass; we considered values of y in the range 





2.0 : 


y = 









Fig. I. — ID velocity dispersion profiles for Dehnen models with- 
out central black holes (solid lines) and with central black holes of mass 
M = [0-^M„^[ (dashed lines), y increases downward, from y = 2 to y = 
(1.5,1.0,0.5)°. 



TABLE 1 
Dehnen model 
parameters 



y 


Op 


0' 


0.5 


0.2775 


0.500 


1.0 


0.3270 


0.330 


1.5 


0.4164 


0.159 


2.0 


0.7071 






0.5 < y < 2. The black hole particle was given zero velocity 
initially. 

Unless otherwise indicated, we present our results in units 
such that G ~ a = Mg^i = 1. In these units, M repre- 
sents the ratio of the black hole mass to the total mass in 
stars. In real galactic spheroids, this ratio is approximately 
10~^, although with so me scatter (Merritt & Ferrarese 20011 
Marconi & Hunri l2003l) . We considered values of M in the 
range 10^^ < M < 10^^ The lower limit is fixed by the re- 
quirement that the black hole particle be appreciably more 
massive than the star particles. 

Figure [l] shows the (ID) radial velocity dispersion profile 
a(r) for Dehnen models with M = and M = 10"^. In the 
absence of a black hole and for y 7^ 2, o(r) peaks at a non-zero 
radius; furthermore the peak value is nearly unaffected by the 
presence of a black hole of mass 10^^Mg„/. Since the am- 
plitude of the Brownian motion of the black hole is expected 
to scale with the stellar velocity dispersion, it is of interest 
to define a characteristic value of a near the centers of these 
models. A natural choice is Op(y), the peak value of o in a 
model without a central black hole. We define rp as the ra- 
dius at which a(r) = a p. Table 1 gives values of Op and rp as 
functions of y. 

Most of our A^-body integrations were carried out on models 
containing stars of equal mass. We also carried out some in- 
tegrations of models in which the stars had a range of masses. 
The stellar masses in the multi-mass integrations were gener- 
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ated randomly from the mass function 

n {m)dm °^ ot^ ' ' ^"^dm (3) 

with a = 1 .35, the Salpeter ( 19551) value. Given a mean mass 
m = Mga]/N with the number of particles, the mass spec- 
trum is determined by equation (|3} and by mi / 1112, the ratio of 
smallest to largest stellar mass. We used mi /m2 = 0.01. 

3. iV-BODY INTEGRATIONS 

All A^-body integrations were carried out using a 
high-accuracy, direct-summation, parallel A^-body code 
((Berczik et al. 2005) on two supercomputers in corporating 
special-purpose GRAPE ( iFukushige et alJ l2005') accelerator 
boards: gravitySimulator' and GRACE.'^ Particle positions 
were advanced using the fourth-order Hermite scheme, with 
discretized individual time steps chosen according to the cri- 
terion of Makino & Aarseth ( 1992), and a time-step accuracy 
parameter r| = 0.01 . All integrations used = 10^ "star" par- 
ticles. Inter-particle forces were softened using the standard 
expression. 



TABLE 2 

Single-mass A'-BODY integrations (A'= 10*) 



'12 ^ 



-Gmim2- 



ri -r2 



|ri-r2| 



3/2 ' 



(4) 



the same softening length was assigned to star and to black 
hole particles. 

In what follows, it is argued that the black hole's Brownian 
motion is determined in large part by perturbations from stars 
that lie within its influence radius r/,, defined (in a slightly 
non-standard way) as the radius containing a mass in stars 
equal to twice the black hole mass. This means that the soft- 
ening length should be small compared with r/,, and in particu- 
lar, small enough that the macroscopic structure of the galaxy 
model does not change at r < r/, due to inaccurate represen- 
tation of the potential. Table 2 gives values of r/, for each set 
of (M,y) values considered here. We verified that a softening 
length of e = 10^^ resulted in no discernible evolution in any 
of the models with y < 1 .5 even at radii <i 0. 1 r/,. Furthermore, 
as shown below (Fig.|3, our results for the y = 1.5 models 
using this value of 8 were essentially indistinguishable from 
those obtained in the earlier study of Dorband et al. (2003) 
who used a code with zero softening. For the models with 
y = 2 and small M, some evolution at r < r/, was observed for 
e = 10^* and we conservatively adopted the smaller softening 
length e = 10^^ for these models. 

Even if e is small enough to accurately reproduce the small- 
scale structure of the galaxy models, evolution can still occur 
after a sufficiently long time as the star particles exchange 
energy. In the equal-mass case, the relevant time scale is the 
two-body relaxation time, 

'"^ ^~ p(r)mG2lnA 

JSpitzedll987i) . with InA die Coulomb logarithm. A colU- 
sional steady state, p oc r~^/^, is reached ne ar the black hole 
T,-{rh) (Bahcall &Woafl976l) . Table 2 gives 
computed using A = r/,a2(r/,)/2Gm* 
(|Preto, Merritt & Spurzem 2004). The short relaxation times 
in models with large y and small M restricts how long these 
models can be integrated before their structure changes. 

' See|http://www.cs.rit.edu/~grapecluster/clus terInfo/grapeCIusterInfo.shtmI| 
^ See^ttpV/www.ari. uni-heidelberg.de/grace 



in a time 
estimates of Triri, 



M 




y 


n, 






Tr(ri,) 


/T/2\ 

{V ) 




Tl 


1 X 10" 


-1 


0.5 


1.11 






6.4 X 10* 


6.72 X 10" 


-6 


2.91 






1.0 


8.09 X 


10" 


-1 


4.3 X 10* 


7.60 X 10" 


-6 


2.37 






1.5 


5.20 X 


10" 


-1 


2.4 X 10* 


1.90 X 10" 


-5 


3.65 






2.0 


2.50 X 


10" 


-1 


9.7 X 10^ 


8.19 X 10" 


-5 


5.46 


3 X 10" 


-2 


0.5 


4.80 X 


; 10" 


-1 


1.4 X 10* 


8.40 X 10" 


-6 


1.09 






1.0 


3.24 X 


10" 


-1 


8.2 X 10' 


2.04 X 10" 


-5 


1.91 






1.5 


1.81 X 


10" 


-1 


3.6 X 10' 


5.74 X 10" 


-5 


3.31 






2.0 


6.38 X 


10" 


-2 


8.4 X 10^ 


2.73 X 10" 


-4 


5.46 


1 X 10" 


-2 


0.5 


2.64 X 


; 10" 


-1 


5.5 X 10' 


2.08 X 10" 


-5 


0.902 






1.0 


1.65 X 


10" 


-1 


2.6 X 10' 


5.49 X 10" 


-5 


1.71 






1.5 


7.95 X 


10" 


-2 


8.2 X 10- 


1.24 X 10" 


-4 


2.39 






2.0 


2.04 X 


10" 


-2 


1.0 X 10- 


7.00 X 10" 


-4 


4.67 


3 X 10" 


-i 


0.5 


1.48 X 


; 10" 


-1 


2.4 X 10^' 


8.87 X 10" 


-5 


1.152 






1.0 


8.40 X 


10" 


-2 


8.3 X 10- 


9.34 X 10" 


-5 


0.873 






1.5 


3.41 X 


10" 


-2 


1.7 X 10- 


2.40 X 10" 


-4 


1.39 






2.0 


6.03 X 


10" 


-3 


1.1 X 10' 


2.13 X 10" 


-3 


4.26 


1 X 10" 


-■i 


0.5 


9.08 X 10" 


-2 


1.2 X 10^ 


1.03 X 10" 


-4 


0.445 






1.0 


4.68 X 


10" 


-2 


3.1 X 10- 


1.96 X 10" 


-4 


0.600 






1.5 


1.61 X 


10" 


-2 


4.1 X 10' 


5.08 X 10" 


-4 


0.977 






2.0 


2.00 X 


10" 


-3 


1.4 X 10° 


5.42 X 10" 


-3 


3.61 


3x 10" 


-4 


0.5 


5.42 X 


; 10" 


-2 


6.6 X 10^ 


1.23 X 10" 


-4 


0.159 






1.0 


2.51 X 


10" 


-2 


1.1 X 10^ 


3.58 X 10" 


-4 


0.335 






1.5 


7.16 X 


10" 


-3 


8.5 X 10" 


1.22 X 10" 


-3 


0.705 






2.0 


6.00 X 


10" 


-4 


1.6 X 10-' 









We adopted the following integration times in model units: 



Y = 


0.5, T 


= 8 


y= 


1.0, T 


= 4 


Y = 


1.5, T 


= 2 


Y = 


2.0, T 


= 1 



As Table 2 shows, the model with (y,M) = (2.0,3 x lO""*) 
has Ti-{rh) ~ 0.16, substantially shorter than the adopted in- 
tegration time, and so short that it would be difficult to get 
good "statistics" on the black hole's Brownian motion in a 
time T,-. Accordingly, we omit the results from this inte- 
gration in the Table and in the discussion that follows. 

When evaluating the amplitude of the black hole's mo- 
tion, a possible concern is any steady component of its ve- 
locity due to a net drift of the A^-body system's center of 
mass. The amplitude of the drift is expected to be of order 
Vmh-VN-^/^, where y is a characteristic internal (stellar) 
velocity. This may be compared with the expected rms veloc- 
ity characterizing the black hole's Brownian motion, which 
is ~ aA'^"'/2(M/Mgai)"'/2. The latter is larger than the for- 
mer by ~ an order of magnitude even for the largest black 
hole mass that we considered, M/Mgai =0.1. We inspected 
the configuration-space motion of the black hole particles in 
all of our integrations; in few if any of the cases could we 
see evidence of a steady drift component to the motion, and 
when it was present, it was too small to significantly affect the 
estimates of (V^). Hence we ignore Vdrfft in what follows. 

The situation is a little different in the multi-mass runs, 
since a single massive particle of mass m^, can spiral into the 
center in a time of order ^ {m/m-^)Ty where m is the mean 
stellar mass. The particle masses in these models were gener- 
ated randomly and a massive particle could find itself initially 
on an orbit for which the dynamical friction time is short. 
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While this effect is a potential source of bias in the multi-mass 
runs, the multi-mass integrations were shorter than the single- 
mass integrations, and we found no evidence of a significant 
time dependence in the Lagrange radii, nor in the black hole's 
(V"), that would suggest a systematic change in the stellar 
distribution over the course of the integrations. In any case, 
the run-to-run variation in (V^) due to different random real- 
izations of the same model were so large that they probably 
swamped this effect (§6). 

4. EXPECTATIONS FROM LOCAL THEORY 

Existing theories of gravitational Brownian motion are too 
idealized to make clear predictions about the amplitude of the 
Brownian motion in our models, for several reasons. First, the 
black hole is massive enough that it modifies the gravitational 
potential as it moves. Second, particles in tightly-bound orbits 
around the black hole will increase its effective mass. Third, 
the properties of the stellar background, e.g. density and ve- 
locity dispersion, are strong functions of distance from the 
black hole. Standard encounter theory (e.g. Spitzer 1987) ap- 
proximates the stellar background as homogeneous and time- 
independent and ignores the effects of a massive particle on 
its surroundings. 

The predictions of idealized theory are nevertheless useful 
as a baseline against which to compare the A^-body results. 
Using the Fokker-Planck equation, the steady-state velocity 
distribution of a massive particle that moves in response to 
perturbations from an infinite and homogeneous distribution 
of background stars is 



fiV)^foe 



(y2) = !^ 



(6) 



where (V^) is the mean square (3D) velocity of the black hole 
and A and C characterize the low-velocity limit of the black 
hole's diffusion coefficients due to encounters with stars: 



(Av|| ) = -Ay+By^ 

{Avj)=C + DV^ ... 



(7a) 
(7b) 



(e.g. Merritt 2001). The latter depend on the phase-space 
number density of stars /j, = /j,(v,r). If /y, is assumed to be 
isotropic in velocity space, then 
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A = —n^G^Mmn 



Jo V 1 



2 ,,4^ 



, (8a) 



C = ^-^n^GWn dvvMv) In (^1 + (8b) 



jMerritll 120051) . Here m is the stellar mass, n is the stel- 
lar number density, and p,„ax is the maximum effective im- 
pact parameter in Chandraskehar's theory. In the special case 
of a Maxwellian distribution of field star velocities, f-k{v) = 

/oe-^'/^"" and 



Ay/2% CMmn , , 

^ = —. ^^W' c = 



8^271 G^m^n 



F{R)^ 



2 Jo 



dze-Hn(l+4RV) 



R=E^^, (9b) 



GM 



and 



^ ' 2A 



(10) 



The predicted (V^) is independent of pmax in this case and is 
equal to the "equipartition" value. 



The velocity distribution at the center of the Dehnen models 
is not precisely Maxwellian. While it is straightforwa rd to de- 
rive expressions for A and C for non-Maxwellian /*'s (iMerrittI 
2005), a more serious problem then presents itself. The veloc- 
ity dispersion at the center of a Dehnen model is zero in the 
absence of a black hole (in all models excepting those with 
y = or 2); and it is infinite when the effect of the black hole 
on the stellar motions is included (Figure[n. Thus equations 
like jlOI I. naively applied, predict either (y^ = or (y^) = oo, 
depending on whether or not the black hole's influence on the 
equilibrium stellar model is taken into account, and of course 
neither result is physically reasonable. The local theory fails 
because the Dehnen models (like real galaxies) are inhomo- 
geneous, both in their density and velocity structure, and be- 
cause the presence of the black hole strongly influences even 
the equilibrium stellar velocity distribution in its vicinity. 

Local theory does make one potentially useful prediction. 
If the massive particle's diffusion coefficients are re-derived 
for the case of a spectrum of perturber masses, one can relate 
the predicted (V^) to its value in the case of a single perturber 
mass. Define the mass function such that n{m)dm is the num- 
ber density of stars in mass range m to m + dm, and assume 
that the velocity distribution is the same for all mass groups. 
In the multi-mass case. 



n{in)mdm 



rdv _piyi^hf_ 

Jo v^*^ >l+pl,y/G^M' 



r _l67t2 2 



n{in)m^dm 







dvvf^,{v) log 



Fmax ^ 
G2M2 



(11a) 



(lib) 



The mean square velocity of the black hole in the multi-mass 
case, for a Maxwellian distribution of field-star velocities, is 
therefore predicted to be 

,2A/_, \M, 



(y2 



(12) 



i.e. m/m times its value in the case of a single perturber mass 
m, where the variable m is 



/ n{m)m^dm 
J n{m)mdm 



(13) 



This prediction will be tested below via A^-body integrations. 

Brownian motion results in a time-dependent displacement 
of the black hole from its otherwise central location in a 
galaxy. The A^-body integrations presented here were judged 
to be too short to extract useful information about the ampli- 
tude of this displacement (as opposed to the velocity), and in 
what follows, we focus on changes in velocity. We will return 
to the question of Brownian displacements in a subsequent 
paper. 

5. RESULTS: SINGLE-MASS MODELS 

FigurelJlillustrates the time dependence of the black hole's 
velocity in a representative set of integrations. The "fre- 
quency" of the motion increases with increasing y, as ex- 
pected, since the period of small oscillations about the poten- 
tial center is smaller for larger y. However none of these plots 
exhibits the quasi-harmonic motion characteristic of particles 
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Fig. 2. — Time dependence of black hole velocity in a set of integrations. From top to bottom, the black hole mass is 10"', 10"^ and 10"^^ in units of the 
galaxy mass. 



wandering in a constant-density core (e.g. Miller & Smith 
1992). 

Brownian motion of a particle in a classical gas obeys a 
Maxwell-Boltzmann distribution. As shown above (§4), the 
same is predicted by the Fokker-Planck equation for a mas- 
sive particle moving in response to gravitational perturbations 
in an infinite homogeneous background of stars, as long as the 
dynamical friction coefficient obeys Hooke's law at low ve- 
locities. However the equations derived above embody many 
approximations that are violated in the A^-body models and in 
real galaxies. Hence it is interesting to directly compute the 
time-averaged velocity distribution of the massive particle in 
the A^-body integrations. Figure|31shows A^(y) for integrations 
with M ~ 10^^ and four different values of y. The velocity of 
the black hole particle was sampled at fixed intervals, ranging 
from At = 0.002 for y ^ 0.5 to Af = 0.00003 for y = 2.0. 

The empirical velocity distributions in Figure |3l are com- 
pared with the Maxwell-Boltzmann distribution, 

N{V)dV = 4%V^ (27i:(y2)/3)"^/^exp {-3V^/2{V^)) dV. 

(14) 

The quantity (V^) in equation (I14> was computed directly 
from the time series of stored velocities; it was not adjusted 
to increase the goodness of fit of equation (I14> to the mea- 
sured A^(y). Neverthless, the fit of the measured A^(y)'s to the 
Maxwell-Boltzmann distribution is clearly good, and nearly 
perfect in the case of y = 2; in this model, the black hole par- 
ticle experienced the largest number of velocity updates. We 
verified that the bumps and wiggles in the plots of N{V) for 
y — 0.5 and y = 1 .0 were essentially random, by constructing 
N{V) separately for the first and second halves of the integra- 
tion intervals. 

If we assume that the time-averaged A^(y) is well described 




5x10"^ 0.01 0.015 0.02 0.01 0.02 0.03 0.04 




0.02 0.04 0.06 0.05 0.1 0.15 0.2 



Fig. 3. — Measured distiibution of black hole velocities in four integra- 
tions with M = 10"^. Thin curves show Maxwell-Boltzmann distributions. 



computed as described in the text. 



by the Maxwell-Boltzmann distribution, then the only quan- 
tity required to specify N{V) is (V^). We expect (V ) to be 
close to the "energy equipartition" value at which the time- 
averaged kinetic energy of the black hole equals the mean ki- 
netic energy of nearby stars. We write this as 



,1, mo 



(15) 



where Op is the peak value of the ID stellar velocity disper- 



6 



Merritt et al. 




M 

Fig. 4. — Mean square velocity of the black hole particle as given in 
Table 2 for the various W-body integrations. Lines show the "equipartition" 
relations, n = 1- R^d diamonds are from Dorband et al. (2003). 



TABLE 3 





iogion = 




A 


+ Blogio 


M 


1 


A 


B 


0.5 


0.844 


0.420 


LO 


0.801 


0.345 


L5 


0.934 


0.307 


2.0 


0.853 


0.093 



sion, as defined in §2, and r| is a parameter whose value is 
expected to be close to one. Values of (V^) are given in Ta- 
ble 12 and plotted in Figures |3 and |5] as functions of M and 
y. These values were computed by sampling the stored veloc- 
ities at small fixed intervals. The lines in Figure |3 show the 
"equipartition" relations, r| = 1 . 
Two results are apparent from Figure s|3 and [S] 

1. There is a clear dependence of (y^) on M, although less 
steep than the expected (V^) °^ M^' dependence. 

2. For a fixed black hole mass, the fractional discrepancy 
between the measured (V^) and the "equipartition" value is 
an increasing function of y. 

Another way to present these results is to interpret equation 
( list as defining the parameter r|, which then measures depar- 
tures from "equipartition" as a function of M and y. Figure |6l 
plots r|, so defined, for these integrations; the r| values are also 
given in Table 2. As the two previous figures indicated, r| in- 
creases systematically with M and y, reaching values as large 
as ^ 5.5 for (M,y) = (0.1,2.0). For small (M,y), r] drops 
below one, to values as low as ^ 0.2. 

The dependence of logrj on logM is nearly linear for fixed 
y (Figure |6j. We carried out least-squares fits to find the co- 
efficients of logioTi w A +Z?logioM. The results of the fits 
are given in Table |3] These results imply that - in a galaxy 
with otherwise fixed properties - the mean square velocity of 
the black hole scales as ^ M^**^ in low-density nuclei and 
~ M^"^ in high-density nuclei. Ignoring the effect of the 
black hole on the stellar motions would lead to the prediction 
(y2) ocM-^ 

There is a natural way to understand these results. The 



steady-state velocities of the stars in these galaxy models are 
influenced by the presence of the black hole (Figure [!), very 
strongly at distances <^ rij. The Brownian motion of the 
black hole depends in turn on perturbations from these fast- 
moving stars. The more massive the black hole, the hotter the 
surrounding stellar fluid, and the greater the expected devi- 
ation of the black hole's Brownian velocity from the value 
that would have been predicted based on the properties of 
a model with the same density structure but lacking a black 
hole. These deviations should also be an increasing function 
of y, since for larger y, the number of high-velocity stars near 
the black hole is larger 

Suppose we assume that the black hole's kinetic energy is 
in equipartition with the stars in some region around it. Define 
to be the ID, mean square stellar velocity within this re- 
gion, which has radius r. The black hole's Brownian velocity 
would then be 



It is understood that in equation ( I16t includes the effect of 
the black hole on the equilibrium stellar motions. We might 
guess, based on the A^-body results, that the best choice for 
f is some multiple of the black hole's radius of influence r/,, 
r = F X ri,. 

Before testing this hypothesis, we make two remarks. 

1. The standard definition of r/j is GM/a^. This definition 
is adequate for galaxies where o(r) is approximately constant 
near the center. For the Dehnen models, this is only true for 
y = 2, and it is never true in models containing a central black 
hole. In §3, we generalized the definition of r/, to be the radius 
containing a mass in stars equal to twice the black hole mass. 
Thus, at r = r/,, the gravitational force from the black hole is 
one-half that due to the stars. We retain that definition in what 
follows. 

2. Above we defined as the mean value of within some 
radius, including the effect of the black hole's gravity on the 
stellar motions. However for y > 2, this quantity diverges as 
a function of the lower integration limit, due to the o= ' 
increase in velocities near the black hole. Nevertheless the 
divergence is only logarithmic as a function of the lower in- 
tegration limit when y = 2. In our A^-body models, a natural 
choice for the lower integration limit is e, the softening length. 
We made this choice in what follows. 

Figure |6l shows the predictions of equation ( I16t compared 
with the A^-body data, setting F = 0.6, i.e. when is com- 
puted as an average within 0.6r/,. In this plot, the ordinate is 
the value of r| that would have been measured if the black 
hole's Brownian velocity were given by equation il6i . i.e. 
r| — a^/a^. The agreement with the data is quite good, con- 
sidering the simplicity of the prescription, the wide range in 
density structure of the models, and the fairly arbitrary defini- 
tion of r/,. The worst fit is for y — 0.5, however these points 
exhibit a large scatter, probably reflecting that the A^-body in- 
tegration time was barely long enough to provide a robust 
sampling of the massive particle's motion (Fig.|2j- 

Other definitions of f were tried (e.g., a fixed value in model 
units) but none was found that reproduced the measured val- 
ues of r| so well. 

We carried out least-squares fits against the data of 
Figure |6l to determine the best value of F for each 
y. For J — (0.5,1,1.5,2), we found best-fit F"s of 
(0.77,0.76,0.52,0.66). The fit to the y = 2 data was nearly 
perfect when the optimal F was used. We conclude that 
the "black holes" in our simulations are in approximate en- 
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ergy equipartition with the stars that lie in a sphere of radius 
(0.67 ± 0. 10)r/, around them. Excluding the rather poor fit for 
Y= 0.5, this becomes (0.65 ±0.10)r/,. 

We can develop a simple model that reproduces these re- 
sults. We suppose that the black hole's motion can be broken 
into two pieces: 

(y2)^(y2) + (y2). (17) 

The first piece, (V^), refers to the motion of the center of mass 
of the system consisting of the black hole and the stars bound 
to it. The second piece, (V^), is the mean square velocity of 
the black hole with respect to the center of mass of the bound 
system. 

Our model requires that we first identify which stars are 
"bound" to the black hole. Clearly, stars with apocenter dis- 
tances r+ <C ri, are bound to the black hole, and stars with 
r+ ^ are not. However stars with r+ w r/, can neither be 
said to be bound or unbound. We will call a star "bound" if 

r+<F'rh. (18) 
Here F' « 1 is a free parameter similar to the factor F defined 
above. According to this definition, the gravitational force on 
the "bound" stars comes predominantly from the black hole. 

Let the mass in bound stars be Mf,. Then the effective mass 
of the bound system (black hole plus stars) is 

Meff=M + Mb. (19) 
The remaining (unbound) stars define a core with some char- 
acteristic density and velocity dispersion. Define the latter to 
be a„. We compute Oj, by taking an average over the unbound 
stars within a region of radius 2F'r/,. 

Given these definitions, the contribution to the black hole's 
motion from the unbound stars is given by Chandrasekhar's 
theory (§4): 



(20) 



In other words, we expect the bound system, of mass Mf,ff, 
to act like a single particle in energy equipartition with the 
unbound stars in the core. 

The motion of the black hole with respect to the center of 
mass of the bound system can be computed by setting 

M\2+m ^ v,=0, (21) 

bound 

where V2 and v, are the velocity of the black hole and of a 
bound star with respect to the center of mass of the bound sys- 
tem. Ignoring correlations between the motions of the bound 
stars, we find 

(V|) = 2^o| (22) 

with (3b the velocity dispersion of the bound stars. 
Combining the two pieces gives 

(y2) 3m/ M 



Since Mf, ; 



M \M + Mt 
M, this is roughly 

^ ' M \ 2 



Mb 2 



(23) 



(24) 



For F' 1 , Of, is of the same order or greater than Gi, for most 
(Y,M), hence the quantity in parentheses on the right hand 
side of equation \2A\ is comparable to 0^. Thus our model 
naturally reproduces the dependence of (V^) observed in the 
A^-body simulations. 

Figure Q shows predicted values of (V^) for F' — \. The 
match with the data in Figure|4]is quite good. 




Fig. 5. — Mean square velocity of the black hole particle for the W-body 
integrations listed in Table 2. These are the same data as in Figure |4| but 
plotted to show the dependence on y. 




M 

Fig. 6. — The quantity T| (equation 15) that measures departures from 
"equipartition" in the W-body models. Lines are the values of T| predicted 
if the black hole is assumed to be in energy equipartition with stars lying at 
distances < 0.6r;, from it. 




Fig. 7. — Predicted values of (V^) based on the model of §5. Symbols 
show the A'-body data, as in Figure 4. 
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Fig. 8. — Mean square black hole velocity in the multi-mass integrations. 
For each value of M, ten W-body integrations were carried out using different 
seeds for the random number generator for the stellar masses; the results for 
each integration are indicated by — . The o symbols show the predicted values 
based on the theory in §4. 



6. RESULTS: MULTI-MASS MODELS 

Stars in real galaxies do not all have the same mass. We 
carried out an additional set of A^-body integrations based on 
the same galaxy models as before, but gen erating the mass of 
each "star" particle from a lSalpeteJ(ll955h mass function, 



n[m)dma:m '■'^"•'t/m, 



mi < m < m2, (25) 

with a = 1.35. We chose mj/nii = 100; the values of mi 
and m2 are then fixed by the requirement that the mean stellar 
mass equal A^^'. Each of these integrations used M = 10^^ 
and = 10^. Because these galaxy models contain a small 
number of massive stars (those at the high-mass tail of the dis- 
tribution), we expect the statistical variance of the results to be 
greater than in the single-mass runs. Hence we carried out 10 
integrations of each model using different seeds to initialize 
the random number generator for the stellar masses. The (V ^) 
values given in Table0]are averages over these 10 runs. We 
also used shorter integration times than in the single-mass in- 
tegrations in order to reduce the effects of mass segregation: 
T = (2,1,0.5,0.25) for y = (0.5,1,1.5,2) respectively. We 
give in Table 0]a^^/2J,, the variance in the measured (V^) val- 
ues for the 10 runs, and ratio between the mean (y2) and the 
value measured in the single-mass runs for the same (y,M,A^). 
We computed the latter quantities over the same time intervals 
used for the multi-mass integrations. 

As shown in §4, Chandrasekhar's theory allows us to pre- 
dict a relation between the black hole's mean square velocity 
in the multi-mass and single-mass cases. The former is pre- 
dicted to be larger than the latter by a factor m/m, where m 
is defined in equation il3\ . Substituting equation i25\ into 
equation (I13> gives 



1- a 

2- a 



mi 



2-a 
m2\ _ ^ 



l-a 

ail) _ 1 



(26) 



Setting a = 1.35 and m2/mi = 100 yields m = ll.lmi = 
4.10m where m — N^^ is the mean stellar mass. Hence, in 
the multi-mass integrations, Chandrasekhar's theory predicts 



TABLE 4 
Multi-mass W-BODY integrations 





y 






ratio 




0.5 


6.58 X IQ-^ 


2.8 X 10-5 


3.98 




1.0 


1.60 X 10-" 


7.0 X 10-5 


3.22 




1.5 


3.77 X 10-" 


1.4 X 10-" 


4.98 




2.0 


1.99 X 10-3 


6.9 X 10-" 


3.87 



a mean square velocity for the black hole that is roughly four 
times larger than in the single-mass integrations. 

The final column of Table |3 gives the measured values of 
this ratio, and Figure|8lshows the 10 measured values of (V^) 
for each value of y, compared with 4.1 times the measured 
value in the corresponding single-mass run. While the vari- 
ation from run to run is large, the mean increase in (V^) is 
reasonably close to the predicted factor. 

7. DISCUSSION 

Our A^-body integrations show that a black hole at the center 
of a dense stellar system responds to perturbations from pass- 
ing stars by achieving a time-averaged kinetic energy given 
by 



3 --2 

-ma 

2 



(27) 



where is the ID, mean square stellar velocity within a re- 
gion r <, 0.6r/, around the black hole, r/, is the black hole's 
influence radius (defined as the radius containing a mass in 
stars equal to twice the black hole's mass), and m is an ef- 
fective stellar mass defined via equation il3\ : in the case of a 
delta-function mass spectrum, m is equal to the stellar mass. 
Equation ( I27> was found to provide a reasonably accurate de- 
scription of the A^-body integrations in galaxy models with 
power-law nuclei, p oc r^y^ 0.5 < y < 2, and for black hole 
masses in the range 10^^ < M/Mgai < 10^'. A more accurate 
description of the A^-body results is given by equation M5\ 
(with m replaced by m) and Table 3, or simply by the mea- 
sured values in Tables 2. The detailed distribution of black 
hole velocities, A^(V), was found to be essentially indistin- 
guishable from a Maxwell-Boltzmann distribution (Figure|3|l. 

As discussed above (§5), equation \21\ has a simple physi- 
cal interpretation: the presence of the black hole deepens the 
galaxy's central potential and increases the velocities of stars 
in its vicinity, and these high-velocity stars provide in turn the 
dominant perturbations that determine the amplitude of the 
black hole's Brownian velocity. The dependence of (V^) on 
M is therefore less steep than ' , since increasing M also 
increases the velocities of nearby stars. 

If we accept equation \21\ . it is straightforward to derive 
simple analytic expressions for the dependence of (V^) on 
M/M gal and y. For the Dehnen models studied here, the black 
hole's influence radius is given by 



rii 



(2M) 



i/(3-r) 



1-(2M) 1/(3-7) 



(2M) (s-T) , M < 1 



(28) 



in model units (G = Mg^i = a = 1). The quantity a that ap- 
pears in equation ( I27t is defined as a number-weighted mean 
within a distance F x r/, from the black hole, where F « 0.6 
(§5). Here we can make use of the analytic expressions for 
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O (r) given bv lDehnenl ( fT99l and lTremaine"etal](IT994l) . We 

write 



(5'{r)=(5\{r)+M(5i{r) 



(29) 



where (3\ is the contribution to o from the stellar potential, 
and O2 is the additional component due to the presence of the 
black hole. For 1 < y < 2 and at radii <, r;,, we have 

, M 1 



af(r): 



2(Y-1)' 



aiir) 



(1+Y) 



(30) 



Averaging over a sphere of radius Fri„ and expressing the re- 
sult in physical units, we find 



V a 



M 



gal 



(Y-2)/(y-3) 



-{l+y){2-y)F^-' 



-(Y-1)(5-2y)' 



(y2-1)(5-2y)(2-y) 



(31a) 



.(31b) 



(Note the divergence when Y= 2, discussed in §5.) Setting 
F = 0.6 (§5) gives // w 2.3 ± 0.4 for Y in the range 1 .2 < Y < 
1.8. The black hole's mean square velocity is then predicted 
to be, via equation ( I27> . 



(y2 



GM 



gal 



M, 



gal 



M 



l/(3-Y) 



1 ^Y^2. 
(32) 



When translated back into model units, this is (V^) w 
7A^"'M"'/(^"^), which is a tolerable fit to the (V^) values 
plotted in Figures 0]and0 The analytic expressions for (V^) 
when Y is in the range 0.5 < Y < 1 are more complicated but 
still imply (V^) oc M-'/^^-T^) for small M/M^ai- 

The fact that (V^) falls off less steeply than M^^ with in- 
creasing mass is perhaps the first clear indication that gravi- 
tational Brownian motion differs in a significant way from its 
fluid analog, for which (V^) °^ ' . It would be of interest to 
refine the M-dependence. Extending the A^-body integrations 
to longer times would reduce the noise but increase the pos- 
sibility of systematic errors due to coUisional evolution of the 
models. An alternative would be to carry out a large number 
of shorter integrations and average the results, as was done 
here in the multi-mass studies (§6). 

When making predictions about the expected amplitude of 
the Brownian velocity of black holes in real galaxies, an ex- 
pression like equation i32\ is not ideal since it contains the 
terms Mgai and a that depend on the large-radius properties 
of the galaxy. A less model-dependent way to present these 
results is in terms of the central properties of the galaxy. Here 
we return to equation il5\ and to the empirical relation estab- 
lished between r| and M: 



logioTi^A+BlogioM 



(33a) 
(33b) 



where the best-fit coefficients A and B as derived from the 
A^-body integrations are given in Table 3. Equations i33\ give 
(y^) in terms of Op, the peak value of a measured outside 
of the black hole's influence radius (Figure[2 Table 1). This 
quantity is easily accessible via ground-based observations for 
many galaxies. 

As an example, we consider the nucleus of the Milky 
Way. The peak velocity dispersion is Op w 150 km s^^ 



jKentlll99^ . For the stellar density profi le at the center of 
the Milky Way bulge. iGenzel et alJ J2003h find p <^ r"^ with 
Y~ 1.4 ±0.1, consistent with the value Y= 1-5 used here in 
some of the A^-body integration s. Finally, th e black hole mass 
is 3.4 ±0.5 X IO^'Mq (Schodel etal«2()03?). Since the Milky 
Way black hole falls on the tight scaling relations defined 
by th e other secure black hol e masses ( Ferrarese & Merri^ 
120001 iMarconi &Hunlll200l . we assume that M/Mg^i ~ 
1 .25 X 10^^, the mean ratio of black hole mass to bul ge mass 
defined by these galaxies (Merritt & Ferrares^ l2001l) . From 
Table 4 or Figure^ we then have r| w 1 .5 and 



VrmiW 0.17 km s 
o 

X 



-1 



1M,T 



1/2 



M 



3.4 X 10*'M,T 



-1/2 



p 



150 km s- 



(34) 



The characteristic stellar mass m (eauationll3> depends on the 
poorly-known mass function f or stars in the Galactic nucleus 
but m is probably of order (IGenzel et al..2003) . Hence we 
predict Vrms « 0.2 km s^ ' . 

Current limits on the motion of the Milky Way black hole 
are based on proper motion measurements of images of Sagit- 
tarius A* (Vp,,, ^ 8 kms"'; Reid et al. (2003)). Such measure- 
ments require observations over an extended period of time; 
in the case of the Reid et al. (2003) study, the baseline was ~ 7 
yr. The measured velocity is derived from the integrated dis- 
placement. If the black hole's velocity changes significantly 
during this time, the distribution of measured velocities will 
be different than the distribution of instantaneous velocities. 
To evaluate the importance of this effect for proper motion 
measurements in the Galactic center, we used our A^-body re- 
sults to compute displacements of the black hole particle over 
finite times, then computed distributions of the proper motion 
velocities that would be measured from these displacements. 
(We note that the time scale over which the massive particle's 
velocity changes is predicted to be independent of the num- 
ber of perturbers for a given mean perturber density, based 
on the equations in §4.) After scaling to the Galactic center, 
we found almost no dependence of the measured (proper mo- 
tion) velocities on the measurement baseline, for times up to 
~ 10^ yr. We therefore base the following discussion on the 
instantaneous velocity distribution A^(V). 

Given that N{V) is expected to be a Maxwell-Boltzmann 
distribution (equation ll4> . we can use the estimate of Vrms to 
compute the probability that a random velocity measurement 
will exceed any value Vq. We first convert the velocity distri- 
bution of equation MA\ into a distribution over 2D velocities 
Vpin in the plane of the sky. Transforming, we find 



N{Vp„,)dypm^{{v^)/^y 



'Vpme-K"/'(^'^dVp„,. 



(35) 



The probability of measuring Vp„, to be greater than Vo, 
given (y2), is then just exp(-3yoV2(V^)). If (V^) = 
(0.17 kms^')^, the expected value for the Milky Way 
black hole, the probability that Vp,„ will exceed Vq — 
(0.1,0.2,0.3,0.4,0.5)kms-i is (0.60,0.13,9.4 x 10^^2.5 x 
10^^,2.3 X 10^^). In order to have a reasonable chance of 
detecting the black hole's motion, a velocity resolution better 
than 0.3 km s^' will be required. 

Once a bona fide measurement of Vp,„ has been made, the 
black hole's velocity can be converted into an estimate of its 
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mass. Here we make use of Bayes's theorem: 



P{M\Vp„,) = 



P{M)P{Vp,„\M) 
JdM'P{M')P{Vp,„\M')' 



(36) 



where P{a\b)da denotes the probability, given b, of measur- 
ing a in the range a to a + da, and P{M)dM is the prior 
probability that M lies in the range M to M + dM. Here 
"prior" means "given all knowledge prior to the measurement 
of Vpm" A completely uninformed prior would be P{M) — 
constant; another standard choice is P{M) <^ M^'; and of 
cour se one c ould use information about M from other sources 
(e.g. ISchodel et al. (200 3)) to construct even more informed 
priors. 

Adopting the maximally-uninformed prior, and using equa- 
tions (l$3t and ( l35t to write P{Vp,„\M), we find the probabiUty 
distribution of M given a measured Vp,,, ~ Vo to be 



P{M\VQ)dM - 



^0 Me-^'^on^^ldM (37) 



and the probabiUty that M > Mq is 



P{M>Mo\Vq)= 1 



Q_ -MoV„'/2,nr]aj, 



(38) 




Fig . 9 . — Probability distribution for tlie mass of the Milky Way black hole, 
based on measurement of a proper motion velocity (solid curves), or based on 
determination of an upper limit to the proper-motion velocity (dashed curves), 
for four values of Vp,,, in km s^ ' . Tick marks indicate median values. 



(In deriving these equations we have ignored the relatively 
weak dependence of r| on M, Table 4. Including that depen- 
dence is straightforward but results in non-analytic expres- 
sions for P.) We define the "best estimate" ?vf of the black 
hole's mass to be the median value from this distribution. This 



IS 



3.36mr|Op 

y2 • 



(39) 



We propose to call this the "Brownian mass estimator." For 



the Milky Way black hole (Op 
estimated mass is 



5W « 2.8 X lO^Mp 



150 km s-\ ^ w 1.5), the 



Vo 



0.2 km s 



(40) 



The 90% confidence intervals on the mass are given by the 
values of Mo for which P{M > Mo\Vo) (0.05,0.95), or 



0.211svf <M<2.82rW. 



Unlike mass estimators based on the statistics of a large 
sample, the uncertainty associated with Brownian mass es- 
timators is irreducible, and it probably makes more sense 
to interpret a measured Vpm as defining a probability dis- 
tribution for M via equation ( I38> than a most-likely value 
via equation ( I40t . For example, if Vpm for the Milky Way 
black hole were measured to be 0.1 km s^', the prob- 
ability that its mass exceeds (10*^, 10^ 1O^)M0 would be 
(99.0%, 56.4%, 0.00058%). 

Figure |9lplots P(M|Vo) in the case of the Milky Way black 
hole (Op = 150 km s^\ r] = 1.5) for four values of Vo- Also 
shown (by tick marks) are the median values ^^{ . 

If we adopt instead the "informed" prior, P{M) oc ^ , then 



P{M\VQ)dM = _}A_.e-^v^-/^''"^4 , (42a) 



P{M > Mo\Vo)=e-^°^o/2>nnoj 



(42b) 



and the Brownian mass estimator becomes 

l.39mx\al _ ,^f... f m \ f Vo 



1.2 X 10°M, 



© 



IMr. 



0.2 km s 



-1 



-2 



(43) 

The median mass is now smaller due to the prior which disfa- 
vors high masses. 

Finally, we ask what can be learned about the black hole 
mass if only an upper limit on its velocity, Vp,,, < Vi,p, is avail- 
able. We consider the determination of an upper limit on Vpm 
to be equivalent to the statement that any velocity in the inter- 
val < Vpm < Vi,p is equally likely, and that velocities greater 
than Vi,p have zero probability. The probability distribution 
for the mass is then 



P{M\V < V, 



Jo 



•p{M\Vo)dVo 



up) 



fo""dVo 



(44) 



(41) Again adopting the uninformed prior for P{M), we find 




-2e-<l^' 



with 



2my\aj, 
M 



(45) 



(46) 



This function is plotted in Figure reffig:prob for various V,ip, 
using the same values of m and Op given above for the Milky 
Way. The low-M tail of the distribution is similar to that of 
P(M|Vo), but there is a more extended tail at high M corre- 
sponding to the fact that low values of Vpm are (by assump- 
tion) equally as likely as Vi,p. 
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